library(pacman)

p_load(tidyverse, marginaleffects, ggplot2, xtable, 
       psych, knitr)


dat <- readRDS("cleaned_svy_dat.Rds")


# Alphas ---------------------------------

# Study 1 #

# Alpha for DVs
dat |> 
  filter(study == 1) |> 
  select(conf_ine_1_recode, conf_ine_2_recode,
         conf_ine_3_recode, conf_ine_4_recode) |> 
  alpha()

dat |> 
  filter(study == 1) |> 
  summarise(mean = mean(conf_ine, na.rm = T),
            sd = sd(conf_ine, na.rm =T),
            missing = sum(is.na(conf_ine)))

# Study 2 #

# Alpha
dat |> 
  filter(study == 2) |> 
  select(conf_ine_1_recode, conf_ine_2_recode,
         conf_ine_3_recode, conf_ine_4_recode) |> 
  alpha()


dat |> 
  filter(study == 2) |> 
  summarise(mean = mean(conf_ine, na.rm = T),
            sd = sd(conf_ine, na.rm =T),
            missing = sum(is.na(conf_ine)))


# Descriptive Sample Statistics ---------------------------

describer <- function(variable, study) {
  
  tab <-  dat |> 
    filter(study == !!study) |>  
    count({{ variable }}) |>  
    mutate(percent = round((n / sum(n)) * 100, 1)) |> 
    na.omit()
  
  return(tab)
  
}

# Study 1 Stats 
for(var in c("gender", "sel", "vote_intent", "age")) {
  
  describer(variable = !!sym(var), study = 1) |> 
    xtable() |> 
    print(include.rownames = F)
  
}


# Study 2 Stats 
for(var in c("gender", "sel", "vote_intent", "age")) {
  
  describer(variable = !!sym(var), study = 2) |> 
    xtable() |> 
    print(include.rownames = F)
  
}




# Confidence in the INE in the Control Condition ------------------

conf_control <- dat |> 
  mutate(
    # minor recodes 
    study_char = if_else(study == 1, "Study 1", "Study 2"),
    support = factor(
      case_match(vote_intent, 
                 "Claudia" ~ "MORENA", 
                 "Xochitl" ~ "PAN/PRI/PRD",
                 "None/Other" ~ "None/Other"),
      levels = c("MORENA",
                 "PAN/PRI/PRD",
                 "None/Other"
      ))
  ) |> 
  # just the control group
  filter(tx == "control") |>
  # drop missings 
  drop_na(conf_ine, support) |> 
  # group the data 
  group_by(support, study_char) |> 
  # means and 95% CIS
  summarise(mean_conf = mean(conf_ine),
            sd = sd(conf_ine), 
            n = n(), 
            se = sd/sqrt(n), 
            upper = mean_conf + (2*se),
            lower = mean_conf - (2*se))  


# plot the means
ggplot(conf_control, aes(x=study_char, 
                         y = mean_conf,
                         ymin = lower,
                         ymax = upper, 
                         color = support)) +
  geom_pointrange() +
  facet_grid(~ support, switch = "x") +  # Ensure to only use one facet_grid function
  labs(y = "Mean Confidence", 
       #title = "Mean Confidence in INE in the Control Condition",
       x = "") +
  theme_classic() +
  theme(
    panel.spacing = unit(0, "lines"),
    strip.placement = "outside",
    panel.border = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(face = "bold"), 
    legend.position = "none", 
    #plot.title = element_text(hjust = 0.5),  
    axis.text = element_text(color = "black")  
  )



ggsave("figures/confidence_control_by_party.pdf", 
       height = 4, width = 4)



# Study 1 model 
party_means_s1 <- lm(conf_ine ~ vote_intent, data = subset(dat, study == 1 &
                                                             tx_char == "control"))
# Study 2 model
party_means_s2 <- lm(conf_ine ~ vote_intent, data = subset(dat, study == 2 &
                                                             tx_char == "control"))

# Comparisons by Party Support - Study 1 
s1_control_means <- avg_comparisons(party_means_s1, 
                                    variable = list(vote_intent = "pairwise")) |> 
  mutate(study = "Study 1")


# Comparisons by Party Support - Study 2 
s2_control_means <- avg_comparisons(party_means_s2, 
                                    variable = list(vote_intent = "pairwise")) |> 
  mutate(study = "Study 2")

# Bind and Print Comparisons  
rbind(s1_control_means, s2_control_means) |> 
  mutate(
    # Recode contrasts
    contrast = case_match(contrast, 
                          "None/Other - Claudia" ~ "None/Other - MORENA",
                          "Xochitl - Claudia" ~ "PAN/PRI/PRD - MORENA",
                          "Xochitl - None/Other" ~  "PAN/PRI/PRD - None Other")
  ) |> 
  # Grab the columns we want
  select( 
    "Study" = study, 
    "Comparison" = contrast, 
    "Estimate" = estimate,
    "Z-Statistic" = statistic, 
    "P-Value" = p.value
  ) |> 
  # Print to latex 
  xtable(digits = 2) |> 
  print(include.rownames = FALSE, 
        file = "tables/party_comparisons_in_control.tex")




# Opinion About Plan B ----------------------------------------


dat <- dat |>
  mutate(
    
    # heard of plan B
    heard_reform = case_match(heard_reform,
                              "No" ~ 0, 
                              "Sí" ~ 1),
    # approve of plan B
    approve_reform = case_match(op_reform, 
                                "Totalmente de acuerdo" ~ 1,
                                "De acuerdo" ~ 1,
                                "Ni de acuerdo ni en desauerdo" ~ 0,
                                "En desacuerdo" ~ 0,
                                "Totalmente en desacuerdo" ~ 0)
  )


# Percent approval of Plan B by Party Support
dat |>  
  group_by(vote_intent) |> 
  summarise(
    heard_of_plan_b = mean(na.omit(heard_reform)),
    approve_of_plan_b = mean(na.omit(approve_reform))
  ) |> 
  na.omit() |> 
  kable(digits = 2)

